##Load libraries and funtions --------------------

{library(lme4)
library(haven)
library(tidyverse)
library(sjPlot)
library(sjstats)
library(jtools)
library(ggplot2)
library(jmv)
library(margins)
library(sjmisc)
library(stats)
library(broom.mixed) 
library(Hmisc)
  library(stargazer)
  library(cowplot)
  library(srvyr)
  library(ggeffects)
  library("marginaleffects")


options(scipen=999)# non-scientific notation

getwd()
setwd("/Users/adamstefkovics/Documents/SCIENCE/ESS/gender_ESS_full")

##loading data

data <- as.data.frame(read_sav('ESS11.sav'))

int_data <- as.data.frame(read_sav('ESS11I.sav'))

cntry_data <- as.data.frame(openxlsx::read.xlsx('cntry_level.xlsx'))

cntry_data %>%
  summarise(
    mean_gii = mean(gii, na.rm = TRUE),
    sd_gii = sd(gii, na.rm = TRUE),
    range_gii = max(gii, na.rm = TRUE) - min(gii, na.rm = TRUE),
    iqr_gii = IQR(gii, na.rm = TRUE),
    
    mean_seats = mean(seats, na.rm = TRUE),
    sd_seats = sd(seats, na.rm = TRUE),
    range_seats = max(seats, na.rm = TRUE) - min(seats, na.rm = TRUE),
    iqr_seats = IQR(seats, na.rm = TRUE)
  )

data$idno_new <- paste(data$cntry, data$idno, sep = "_")
int_data$idno_new <- paste(int_data$cntry, int_data$idno, sep = "_")
int_data$intnum_new <- paste(int_data$cntry, int_data$intnum, sep = "_")

##joining
data <- data %>% 
  left_join(select(int_data, -cntry,-essround), by = "idno_new")

data <- merge(data,cntry_data,by="cntry",all.x = TRUE)

###exclude certain interviewers
int_num <- data %>% 
  group_by(intnum_new) %>% 
  summarise(int_num=n()) %>% 
  filter(int_num<10)


data <- merge(data,int_num,by="intnum_new",all.x=T) %>% 
  filter(is.na(int_num))

gender_num <- data %>% 
  group_by(intnum_new,gndr) %>% 
  summarise(gender_case=n(),.groups = "drop") %>% 
  complete(intnum_new, gndr, fill = list(gender_case = 0)) %>% 
  filter(gndr==1) %>% 
  select(intnum_new,gender_case)


data <- merge(data,gender_num,by="intnum_new",all.x=T) %>% 
  filter(gender_case>0)

data <- data %>% filter(!is.na(intgndr))

data %>% 
  group_by(intnum_new) %>% 
  summarise(n())

###transform

data <- data %>% 
  mutate(gender_health=case_when(trmdcnt==1~1,
                                 trmdcnt>1~0,
                                 TRUE~NA_real_),
         gender_work=case_when(trwkcnt==1~1,
                               trwkcnt>1~0,
                               TRUE~NA_real_),
         
         gender_police=case_when(trplcnt==1~1,
                                 trplcnt>1~0,
                                 TRUE~NA_real_),
         gender=gndr,
         gender_int=intgndr,
         marital=case_when(rshpsts==1~1,
                           rshpsts==4~1,
                           TRUE~0),
         work=case_when(pdwrk==1~1,
                        TRUE~0),
         child = case_when(
           rshipa2 == 2 | rshipa3 == 2 | rshipa4 == 2 | rshipa5 == 2 | 
             rshipa6 == 2 | rshipa7 == 2 | rshipa8 == 2 | rshipa9 == 2 | 
             rshipa10 == 2 | rshipa11 == 2 | rshipa12 == 2 | chldhhe==1 ~ 1,
           TRUE ~ 0),
         christian=case_when(rlgdnm<5~1,
                             TRUE~0),
         unfair_med=case_when(trmedmw<3~1,
                              TRUE~0),
         unfair_med_inr=case_when(is.na(unfair_med)~0,
                              TRUE~1),
         unfair_work=case_when(trwrkmw<3~1,
                               TRUE~0),
         unfair_work_inr=case_when(is.na(unfair_work)~0,
                              TRUE~1),
         unfair_police=case_when(trplcmw<3~1,
                                 TRUE~0),
         unfair_police_inr=case_when(is.na(unfair_police)~0,
                               TRUE~1),
         domicil_3=as.factor(case_when(domicil<3~3,
                                     domicil==3~2,
                                     domicil>3~1,
                                     TRUE~NA_real_)),
         unfair=case_when(unfair_med==1 ~1,
                                unfair_work==1~1,
                                unfair_police==1~1,
                              TRUE~0),
         relig=rlgdgr,
         edu_cat=case_when(edulvlb<313~1,
                           edulvlb %in% 313:423~2,
                           edulvlb > 423~3,
                           TRUE~NA_real_),
         edulvlb_cont = case_when(
             edulvlb == 0    ~ 1,
             edulvlb == 113  ~ 2,
             edulvlb == 129  ~ 3,
             edulvlb == 212  ~ 4,
             edulvlb == 213  ~ 5,
             edulvlb == 221  ~ 6,
             edulvlb == 222  ~ 7,
             edulvlb == 223  ~ 8,
             edulvlb == 229  ~ 9,
             edulvlb == 311  ~ 10,
             edulvlb == 312  ~ 11,
             edulvlb == 313  ~ 12,
             edulvlb == 321  ~ 13,
             edulvlb == 322  ~ 14,
             edulvlb == 323  ~ 15,
             edulvlb == 412  ~ 16,
             edulvlb == 413  ~ 17,
             edulvlb == 421  ~ 18,
             edulvlb == 422  ~ 19,
             edulvlb == 423  ~ 20,
             edulvlb == 510  ~ 21,
             edulvlb == 520  ~ 22,
             edulvlb == 610  ~ 23,
             edulvlb == 620  ~ 24,
             edulvlb == 710  ~ 25,
             edulvlb == 720  ~ 26,
             edulvlb == 800  ~ 27,
             edulvlb == 5555 ~ NA_real_,  # Handle Other/Missing values as NA
             edulvlb == 7777 ~ NA_real_,
             edulvlb == 8888 ~ NA_real_,
             edulvlb == 9999 ~ NA_real_,
             TRUE ~ NA_real_),  # Default to NA for all unrecognized values
         relig_reverse = case_when(
           rlgatnd == 1 ~ 7,
           rlgatnd == 2 ~ 6,
           rlgatnd == 3 ~ 5,
           rlgatnd == 4 ~ 4,
           rlgatnd == 5 ~ 3,
           rlgatnd == 6 ~ 2,
           rlgatnd == 7 ~ 1,
           TRUE ~ NA_real_),
         hincfel_rev = case_when(
           hincfel == 1 ~ 4,  
           hincfel == 2 ~ 3,  
           hincfel == 3 ~ 2,  
           hincfel == 4 ~ 1,
           TRUE ~ NA_real_),
         prewhp_new=case_when(is.na(prewhp)~0,
                              TRUE~prewhp),
         eqwrkbg_inr=case_when(is.na(prewhp)~0,
                               TRUE~1),
         eqpolbg_inr=case_when(is.na(eqpolbg)~0,
                               TRUE~1),
         eqmgmbg_inr=case_when(is.na(eqmgmbg)~0,
                               TRUE~1),
         eqpaybg_inr=case_when(is.na(eqpaybg)~0,
                               TRUE~1),
         eqparlv_inr=case_when(is.na(eqparlv)~0,
                               TRUE~1),
         eqparep_inr=case_when(is.na(eqparep)~0,
                               TRUE~1),
         freinsw_inr=case_when(is.na(freinsw)~0,
                               TRUE~1),
         fineqpy_inr=case_when(is.na(fineqpy)~0,
                               TRUE~1),
         wsekpwr_inr=case_when(is.na(wsekpwr)~0,
                               TRUE~1),
         weasoff_inr=case_when(is.na(weasoff)~0,
                               TRUE~1),
         wexashr_inr=case_when(is.na(wexashr)~0,
                               TRUE~1),
         eqparlv_rev = case_when(
           eqparlv == 1 ~ 5,  
           eqparlv == 2 ~ 4, 
           eqparlv == 3 ~ 3,  
           eqparlv == 4 ~ 2,  
           eqparlv == 5 ~ 1,  
           eqparlv %in% c(7, 8, 9) ~ NA_real_,  
           TRUE ~ NA_real_),
        eqparep_rev = case_when(
             eqparep == 1 ~ 5,  
             eqparep == 2 ~ 4, 
             eqparep == 3 ~ 3,  
             eqparep == 4 ~ 2,  
             eqparep == 5 ~ 1,  
             eqparep %in% c(7, 8, 9) ~ NA_real_,  
             TRUE ~ NA_real_),
        freinsw_rev = case_when(
          freinsw == 1 ~ 5,  
          freinsw == 2 ~ 4, 
          freinsw == 3 ~ 3,  
          freinsw == 4 ~ 2,  
          freinsw == 5 ~ 1,  
          freinsw %in% c(7, 8, 9) ~ NA_real_,  
          TRUE ~ NA_real_),
        fineqpy_rev = case_when(
          fineqpy == 1 ~ 5,  
          fineqpy == 2 ~ 4, 
          fineqpy == 3 ~ 3,  
          fineqpy == 4 ~ 2,  
          fineqpy == 5 ~ 1,  
          fineqpy %in% c(7, 8, 9) ~ NA_real_,  
          TRUE ~ NA_real_),
        wsekpwr_rev = case_when(
          wsekpwr == 1 ~ 5,  
          wsekpwr == 2 ~ 4, 
          wsekpwr == 3 ~ 3,  
          wsekpwr == 4 ~ 2,  
          wsekpwr == 5 ~ 1,  
          wsekpwr %in% c(7, 8, 9) ~ NA_real_,  
          TRUE ~ NA_real_),
        weasoff_rev = case_when(
          weasoff == 1 ~ 5,  
          weasoff == 2 ~ 4, 
          weasoff == 3 ~ 3,  
          weasoff == 4 ~ 2,  
          weasoff == 5 ~ 1,  
          weasoff %in% c(7, 8, 9) ~ NA_real_,  
          TRUE ~ NA_real_),
        wexashr_rev = case_when(
          wexashr == 1 ~ 5,  
          wexashr == 2 ~ 4, 
          wexashr == 3 ~ 3,  
          wexashr == 4 ~ 2,  
          wexashr == 5 ~ 1,  
          wexashr %in% c(7, 8, 9) ~ NA_real_,  
          TRUE ~ NA_real_),
        vdcond_home=case_when(vdcond==1~1,
                              TRUE~0),
        vdcond_out=case_when(vdcond==2~1,
                              TRUE~0),
        vdcond_video=case_when(vdcond==3~1,
                              TRUE~0))

##function for scaling
scale01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))} 

data$gender <- as.factor(data$gender)
data$gender_int <- as.factor(data$gender_int)
data$vdcond <- as.factor(data$vdcond)
data$prewhp_new <- as.factor(data$prewhp_new)
data$lrscale_sc = scale01(data$lrscale)
data$lrscale_c = scale(data$lrscale,center=TRUE,scale=FALSE)
data$age_sc = scale01(data$agea)
data$age_c = scale(data$agea,center=TRUE,scale=FALSE)
data$intage_sc = scale01(data$intagea)
data$intage_c = scale(data$intagea,center=TRUE,scale=FALSE)
data$relig_reverse_sc = scale01(data$relig_reverse)
data$relig_reverse_c = scale(data$relig_reverse_sc,center=TRUE,scale=FALSE)
data$edulvlb_cont_sc = scale01(data$edulvlb_cont)
data$edulvlb_cont_c = scale(data$edulvlb_cont,center=TRUE,scale=FALSE)
data$gggi_sc = scale01(data$gggi)
data$gii_sc = scale01(data$gii)
data$mortality_sc = scale01(data$mortality)
data$birth_sc = scale01(data$birth)
data$seats_sc = scale01(data$seats)
data$labor_sc = scale01(data$labor)
data$vdcond_home <- as.factor(data$vdcond_home)
data$vdcond_out <- as.factor(data$vdcond_out)
data$vdcond_video <- as.factor(data$vdcond_video)
data$vdcond <- as.factor(data$vdcond)
data$unfair_med<-as.factor(data$unfair_med) 
data$unfair_work<-as.factor(data$unfair_work) 
data$unfair_police<-as.factor(data$unfair_police) 
data$unfair<-as.factor(data$unfair) 


##creating gender dyads
data <- data %>% 
  mutate(dyad=case_when(gender==1 & gender_int==1~1,
                        gender==1 & gender_int==2~2,
                        gender==2 & gender_int==1~3,
                        gender==2 & gender_int==2~4))

frq(data,dyad)
data$dyad <- as.factor(data$dyad)


##exclude where PSU is not uniqe

country_duplicate_flag <- data %>%
  group_by(cntry) %>%
  summarise(has_duplicates = ifelse(n() > n_distinct(psu), 1, 0)) 

# Merge the flag back to the original dataset
data <- data %>%
  left_join(country_duplicate_flag, by = "cntry") %>% 
  filter(has_duplicates == 1)


}


##analyis


###H1a

dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")

coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the logistic regression formula dynamically
    formula <- as.formula(paste(dv, "~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                child + work + hincfel + relig_reverse_sc + lrscale_sc + gender_int + intage_sc +
                vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry)"))
    
    # Fit the logistic mixed-effects model
    model <- glmer(formula, data = data, family = binomial, weights = pspwght,
                   control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
    
    # Extract coefficients, standard errors, and z-values
    coef_summary <- summary(model)$coefficients
    dv_int_row <- grep(paste0("^", "gender_int2", "$"), rownames(coef_summary))
    
    gender_int_coefficient <- coef_summary[dv_int_row, "Estimate"]
    gender_int_std_error <- coef_summary[dv_int_row, "Std. Error"]
    gender_int_z_value <- coef_summary[dv_int_row, "z value"]
    
    # Calculate the p-value
    gender_int_p_value <- 2 * (1 - pnorm(abs(gender_int_z_value)))
    
    # Create a dataframe with the dependent variable, coefficient, standard error, and p-value
    coef_df <- data.frame(dependent_variable = dv, 
                          Coefficient = gender_int_coefficient, 
                          Std_Error = gender_int_std_error, 
                          P_Value = gender_int_p_value)
    
    # Store the data frame in the list
    coef_list[[dv]] <- coef_df
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}

results_h1a <- coeff_extract(dependent_vars)

openxlsx::write.xlsx(results_h1a,"h1a.xlsx")


##plot H1a
labels=c("Medical treatment",
         "Work",
         "Police encountrer",
         "Any domain")

a <-results_h1a %>% 
  mutate(p_value = case_when(P_Value < 0.001 ~ "1",TRUE ~ "0")) %>% 
  mutate(dependent_variable = factor(dependent_variable, 
                                     levels = c("unfair_med",
                                                "unfair_work",
                                                "unfair_police",
                                                "unfair"))) %>% 
  ggplot(aes(x = dependent_variable, y = Coefficient)) +
  geom_pointrange(aes(ymin = Coefficient - 2 * Std_Error, 
                      ymax = Coefficient + 2 * Std_Error, 
                      color = as.factor(p_value)),
                  linewidth=1.1) +  # Removed `fill` from here
  geom_point(aes(fill = as.factor(p_value),color = as.factor(p_value)), shape = 21, size = 3) +  # Fill only in points
  geom_hline(yintercept = 0, color = "grey75", linetype = "dashed") +
  scale_fill_manual(values = c("grey75", "#688cdb")) +
  scale_color_manual(values = c("grey75", "#688cdb")) +
  theme_minimal() +
  coord_flip() +
  ggtitle("Disclosure")+
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(name = "Est. of the interviewer's gender\non disclosure") +
  scale_x_discrete(limits=rev,
                   expand = expansion(add = 1.2),
                   name = "",
                   labels = rev(labels))

#H1b

dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")


coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the logistic regression formula dynamically
    formula <- as.formula(paste(dv, "~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                child + work + hincfel + relig_reverse_sc + lrscale_sc + 
               gender_int * gender + intage_sc +
                vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry)"))
    
    # Fit the logistic mixed-effects model
    model <- tryCatch(
      glmer(formula, data = data, family = binomial, weights = pspwght,
            control = glmerControl(optimizer = "bobyqa", 
                                   optCtrl = list(maxfun = 100000))),
      error = function(e) {
        message(paste("Error in model fitting for:", dv))
        return(NULL)
      }
    )
    
    # Skip if model fitting failed
    if (is.null(model)) next
    
    # Extract coefficients table
    coef_summary <- summary(model)$coefficients
    
    # Debugging: Print available coefficient names
    print(rownames(coef_summary))
    
    # Adjust grep to be more flexible
    dyad_coef_rows <- grep("gender2:gender_int2", rownames(coef_summary))
    
    # Check if coefficients exist before extracting
    if (length(dyad_coef_rows) > 0) {
      dyad_coefs <- coef_summary[dyad_coef_rows, "Estimate"]
      dyad_std_errors <- coef_summary[dyad_coef_rows, "Std. Error"]
      dyad_z_values <- coef_summary[dyad_coef_rows, "z value"]
      
      # Calculate p-values
      dyad_p_values <- 2 * (1 - pnorm(abs(dyad_z_values)))
      
      # Create a dataframe with the extracted coefficients
      coef_df <- data.frame(
        dependent_variable = rep(dv, length(dyad_coefs)), 
        dyad_category = rownames(coef_summary)[dyad_coef_rows], 
        Coefficient = dyad_coefs, 
        Std_Error = dyad_std_errors, 
        P_Value = dyad_p_values
      )
      
      # Store the data frame in the list
      coef_list[[dv]] <- coef_df
    } else {
      message("No dyad coefficients found for: ", dv)
    }
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}

results_h1b <- coeff_extract(dependent_vars)

results_h1b
openxlsx::write.xlsx(results_h1b,"h1b.xlsx")


##H2a

dependent_vars <- c("eqwrkbg",
                           "eqpolbg",
                           "eqmgmbg",
                           "eqpaybg",
                           "eqparep_rev",
                           "eqparlv_rev",
                           "freinsw_rev",
                           "fineqpy_rev",
                           "wsekpwr_rev",
                           "weasoff_rev",
                           "wexashr_rev")

coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the formula dynamically, using dv_int as the independent variable
    formula <- as.formula(paste(dv, "~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry)"))
    
    # Fit the model
    model <- lmer(formula, data = data, REML = FALSE,weights=pspwght, 
                  control = lmerControl(optimizer = "optimx", 
                                        optCtrl = list(method = 'nlminb')))
    
    # Extract coefficients, standard errors, and t-values
    coef_summary <- summary(model)$coefficients
    dv_int_row <- grep(paste0("^", "gender_int2", "$"), rownames(coef_summary))
    
    gender_int_coefficient <- coef_summary[dv_int_row, "Estimate"]
    gender_int_std_error <- coef_summary[dv_int_row, "Std. Error"]
    gender_int_t_value <- coef_summary[dv_int_row, "t value"]
    
    # Calculate the degrees of freedom
    df <- nobs(model) - length(fixef(model))
    
    # Calculate the p-value
    gender_int_p_value <- 2 * (1 - pt(abs(gender_int_t_value), df))
    
    # Create a dataframe with the dependent variable, its coefficient, standard error, and p-value
    coef_df <- data.frame(dependent_variable = dv, 
                          Coefficient = gender_int_coefficient, 
                          Std_Error = gender_int_std_error, 
                          P_Value = gender_int_p_value)
    
    # Store the data frame in the list
    coef_list[[dv]] <- coef_df
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}


results_h2a <- coeff_extract(dependent_vars)

labels=c("Equality - family life",
"Equality - politics",
"Equality - work",
"Equality - economy",
"Policy - parliament",
"Policy - parental leave",
"Policy - insult",
"Policy - payments",
"Stereotype - control",
"Stereotype - offend",
"Stereotype - exaggerate")

b<- results_h2a %>% 
mutate(p_value = case_when(P_Value < 0.05 ~ "1",TRUE ~ "0")) %>% 
  mutate(dependent_variable = factor(dependent_variable, levels = c("eqwrkbg",
                                            "eqpolbg",
                                            "eqmgmbg",
                                            "eqpaybg",
                                            "eqparep_rev",
                                            "eqparlv_rev",
                                            "freinsw_rev",
                                            "fineqpy_rev",
                                            "wsekpwr_rev",
                                            "weasoff_rev",
                                            "wexashr_rev"))) %>% 
  ggplot(aes(x = dependent_variable, y = Coefficient)) +
  geom_pointrange(aes(ymin = Coefficient - 2 * Std_Error, 
                      ymax = Coefficient + 2 * Std_Error, 
                  color = as.factor(p_value)),
                  linewidth=1.1) +  # Removed `fill` from here
  geom_point(aes(fill = as.factor(p_value),color = as.factor(p_value)), shape = 21, size = 3) +  # Fill only in points
  geom_hline(yintercept = 0, color = "grey75", linetype = "dashed") +
  scale_fill_manual(values = c("grey75", "#688cdb")) +
  scale_color_manual(values = c("grey75", "#688cdb")) +
    theme_minimal() +
  coord_flip() +
  ggtitle("Distributions")+
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(name = "Est. of the interviewer's gender\non gender attitudes") +
  scale_x_discrete(limits=rev,
                   expand = expansion(add = 1.2),
                   name = "",
                   labels = rev(labels))

p<- plot_grid(a,b,nrow=1,labels=c("a","b"))

p
ggsave("fig1.png",p,h=5,w=8)

##H2b##
  
dependent_vars <- c("eqwrkbg",
                      "eqpolbg",
                      "eqmgmbg",
                      "eqpaybg",
                      "eqparep_rev",
                      "eqparlv_rev",
                      "freinsw_rev",
                      "fineqpy_rev",
                      "wsekpwr_rev",
                      "weasoff_rev",
                      "wexashr_rev")
  
coeff_extract <- function(dependent_vars) {
    # Initialize a list to store coefficient data frames for each dependent variable
    coef_list <- list()
    
    # Loop through each dependent variable
    for (dv in dependent_vars) {
      
      # Construct the formula dynamically, using dv_int as the independent variable
      formula <- as.formula(paste(dv, "~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                child+work+hincfel+relig_reverse_sc+lrscale_sc+
                gender_int*gender+intage_sc+
                vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry)"))
      
      # Fit the model
      model <- lmer(formula, data = data, REML = FALSE,weights=pspwght, 
                    control = lmerControl(optimizer = "optimx", 
                                          optCtrl = list(method = 'nlminb')))
      
      # Extract coefficients, standard errors, and t-values
      coef_summary <- summary(model)$coefficients
      dv_int_row <- grep(paste0("^", "gender2:gender_int2", "$"), rownames(coef_summary))
      
      gender_int_coefficient <- coef_summary[dv_int_row, "Estimate"]
      gender_int_std_error <- coef_summary[dv_int_row, "Std. Error"]
      gender_int_t_value <- coef_summary[dv_int_row, "t value"]
      
      # Calculate the degrees of freedom
      df <- nobs(model) - length(fixef(model))
      
      # Calculate the p-value
      gender_int_p_value <- 2 * (1 - pt(abs(gender_int_t_value), df))
      
      # Create a dataframe with the dependent variable, its coefficient, standard error, and p-value
      coef_df <- data.frame(dependent_variable = dv, 
                            Coefficient = gender_int_coefficient, 
                            Std_Error = gender_int_std_error, 
                            P_Value = gender_int_p_value)
      
      # Store the data frame in the list
      coef_list[[dv]] <- coef_df
    }
    
    # Combine all coefficient data frames into one
    combined_coef_df <- do.call(rbind, coef_list)
    
    return(combined_coef_df)
  }
  
  
results_h2b <- coeff_extract(dependent_vars)

results_h2b  

####################
##H3
####################

dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")


coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the logistic regression formula dynamically
    formula <- as.formula(paste(dv, "~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                child + work + hincfel + relig_reverse_sc + lrscale_sc + 
               gender_int*gii_sc + intage_sc +
                vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry)"))
    
    # Fit the logistic mixed-effects model
    model <- tryCatch(
      glmer(formula, data = data, family = binomial, weights = pspwght,
            control = glmerControl(optimizer = "bobyqa", 
                                   optCtrl = list(maxfun = 100000))),
      error = function(e) {
        message(paste("Error in model fitting for:", dv))
        return(NULL)
      }
    )
    
    # Skip if model fitting failed
    if (is.null(model)) next
    
    # Extract coefficients table
    coef_summary <- summary(model)$coefficients
    
    # Debugging: Print available coefficient names
    print(rownames(coef_summary))
    
    # Adjust grep to be more flexible
    dyad_coef_rows <- grep("gender_int2:gii_sc", rownames(coef_summary))
    
    # Check if coefficients exist before extracting
    if (length(dyad_coef_rows) > 0) {
      dyad_coefs <- coef_summary[dyad_coef_rows, "Estimate"]
      dyad_std_errors <- coef_summary[dyad_coef_rows, "Std. Error"]
      dyad_z_values <- coef_summary[dyad_coef_rows, "z value"]
      
      # Calculate p-values
      dyad_p_values <- 2 * (1 - pnorm(abs(dyad_z_values)))
      
      # Create a dataframe with the extracted coefficients
      coef_df <- data.frame(
        dependent_variable = rep(dv, length(dyad_coefs)), 
        dyad_category = rownames(coef_summary)[dyad_coef_rows], 
        Coefficient = dyad_coefs, 
        Std_Error = dyad_std_errors, 
        P_Value = dyad_p_values
      )
      
      # Store the data frame in the list
      coef_list[[dv]] <- coef_df
    } else {
      message("No dyad coefficients found for: ", dv)
    }
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}

results_h3a <- coeff_extract(dependent_vars)
results_h3a
openxlsx::write.xlsx(results_h3a,"h3a.xlsx")


coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the logistic regression formula dynamically
    formula <- as.formula(paste(dv, "~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                child + work + hincfel + relig_reverse_sc + lrscale_sc + 
               gender_int*seats_sc + intage_sc +
                vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry)"))
    
    # Fit the logistic mixed-effects model
    model <- tryCatch(
      glmer(formula, data = data, family = binomial, weights = pspwght,
            control = glmerControl(optimizer = "bobyqa", 
                                   optCtrl = list(maxfun = 100000))),
      error = function(e) {
        message(paste("Error in model fitting for:", dv))
        return(NULL)
      }
    )
    
    # Skip if model fitting failed
    if (is.null(model)) next
    
    # Extract coefficients table
    coef_summary <- summary(model)$coefficients
    
    # Debugging: Print available coefficient names
    print(rownames(coef_summary))
    
    # Adjust grep to be more flexible
    dyad_coef_rows <- grep("gender_int2:seats_sc", rownames(coef_summary))
    
    # Check if coefficients exist before extracting
    if (length(dyad_coef_rows) > 0) {
      dyad_coefs <- coef_summary[dyad_coef_rows, "Estimate"]
      dyad_std_errors <- coef_summary[dyad_coef_rows, "Std. Error"]
      dyad_z_values <- coef_summary[dyad_coef_rows, "z value"]
      
      # Calculate p-values
      dyad_p_values <- 2 * (1 - pnorm(abs(dyad_z_values)))
      
      # Create a dataframe with the extracted coefficients
      coef_df <- data.frame(
        dependent_variable = rep(dv, length(dyad_coefs)), 
        dyad_category = rownames(coef_summary)[dyad_coef_rows], 
        Coefficient = dyad_coefs, 
        Std_Error = dyad_std_errors, 
        P_Value = dyad_p_values
      )
      
      # Store the data frame in the list
      coef_list[[dv]] <- coef_df
    } else {
      message("No dyad coefficients found for: ", dv)
    }
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}

results_h3ab <- coeff_extract(dependent_vars)
results_h3ab
openxlsx::write.xlsx(results_h3ab,"h3ab.xlsx")



##H3b##

dependent_vars <- c("eqwrkbg",
                    "eqpolbg",
                    "eqmgmbg",
                    "eqpaybg",
                    "eqparep_rev",
                    "eqparlv_rev",
                    "freinsw_rev",
                    "fineqpy_rev",
                    "wsekpwr_rev",
                    "weasoff_rev",
                    "wexashr_rev")

coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the formula dynamically, using dv_int as the independent variable
    formula <- as.formula(paste(dv, "~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                child+work+hincfel+relig_reverse_sc+lrscale_sc+
                gender_int*gii_sc+intage_sc+
                vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry)"))
    
    # Fit the model
    model <- lmer(formula, data = data, REML = FALSE,weights=pspwght, 
                  control = lmerControl(optimizer = "optimx", 
                                        optCtrl = list(method = 'nlminb')))
    
    # Extract coefficients, standard errors, and t-values
    coef_summary <- summary(model)$coefficients
    dv_int_row <- grep(paste0("^", "gender_int2:gii_sc", "$"), rownames(coef_summary))
    
    gender_int_coefficient <- coef_summary[dv_int_row, "Estimate"]
    gender_int_std_error <- coef_summary[dv_int_row, "Std. Error"]
    gender_int_t_value <- coef_summary[dv_int_row, "t value"]
    
    # Calculate the degrees of freedom
    df <- nobs(model) - length(fixef(model))
    
    # Calculate the p-value
    gender_int_p_value <- 2 * (1 - pt(abs(gender_int_t_value), df))
    
    # Create a dataframe with the dependent variable, its coefficient, standard error, and p-value
    coef_df <- data.frame(dependent_variable = dv, 
                          Coefficient = gender_int_coefficient, 
                          Std_Error = gender_int_std_error, 
                          P_Value = gender_int_p_value)
    
    # Store the data frame in the list
    coef_list[[dv]] <- coef_df
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}


results_h3b <- coeff_extract(dependent_vars)
results_h3b


coeff_extract <- function(dependent_vars) {
  # Initialize a list to store coefficient data frames for each dependent variable
  coef_list <- list()
  
  # Loop through each dependent variable
  for (dv in dependent_vars) {
    
    # Construct the formula dynamically, using dv_int as the independent variable
    formula <- as.formula(paste(dv, "~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                child+work+hincfel+relig_reverse_sc+lrscale_sc+
                gender_int*seats_sc+intage_sc+
                vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry)"))
    
    # Fit the model
    model <- lmer(formula, data = data, REML = FALSE,weights=pspwght, 
                  control = lmerControl(optimizer = "optimx", 
                                        optCtrl = list(method = 'nlminb')))
    
    # Extract coefficients, standard errors, and t-values
    coef_summary <- summary(model)$coefficients
    dv_int_row <- grep(paste0("^", "gender_int2:seats_sc", "$"), rownames(coef_summary))
    
    gender_int_coefficient <- coef_summary[dv_int_row, "Estimate"]
    gender_int_std_error <- coef_summary[dv_int_row, "Std. Error"]
    gender_int_t_value <- coef_summary[dv_int_row, "t value"]
    
    # Calculate the degrees of freedom
    df <- nobs(model) - length(fixef(model))
    
    # Calculate the p-value
    gender_int_p_value <- 2 * (1 - pt(abs(gender_int_t_value), df))
    
    # Create a dataframe with the dependent variable, its coefficient, standard error, and p-value
    coef_df <- data.frame(dependent_variable = dv, 
                          Coefficient = gender_int_coefficient, 
                          Std_Error = gender_int_std_error, 
                          P_Value = gender_int_p_value)
    
    # Store the data frame in the list
    coef_list[[dv]] <- coef_df
  }
  
  # Combine all coefficient data frames into one
  combined_coef_df <- do.call(rbind, coef_list)
  
  return(combined_coef_df)
}


results_h3c <- coeff_extract(dependent_vars)



###extracting tables

##H1a

dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")

reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence")

model_1<- glmer(unfair_med~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                child + work + hincfel + relig_reverse_sc + lrscale_sc + gender_int + 
                  intage_sc +
                vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry),
      data = data, family = binomial, weights = pspwght,
      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_2<- glmer(unfair_work~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + gender_int + intage_sc +
                  vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_3<- glmer(unfair_police~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + gender_int + intage_sc +
                  vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_4<- glmer(unfair~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + gender_int + intage_sc +
                  vdcond + prewhp_new + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

tab_model(model_1, model_2,model_3,model_4,pred.labels=reg_label,dv.labels=c("Medical treatment",
                                            "Work",
                                            "Police encounter",
                                            "Any domain"),
          collapse.ci=T,file="h1a.doc")


##H1b

dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")

reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence","R's gender X I's gender")

model_1<- glmer(unfair_med~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gender + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

sub_data <- model.frame(model_1)
sub_data <- as.data.frame(model.frame(model_1))
flat_table(sub_data, gender, gender_int)

model_2<- glmer(unfair_work~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gender + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_3<- glmer(unfair_police~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gender+ (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_4<- glmer(unfair~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gender + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Medical treatment",
                                                                          "Work",
                                                                          "Police encounter",
                                                                          "Any domain"),
          collapse.ci=T,file="h1b.doc")

##H2a

dependent_vars <- c("eqwrkbg",
                    "eqpolbg",
                    "eqmgmbg",
                    "eqpaybg",
                    "eqparep_rev",
                    "eqparlv_rev",
                    "freinsw_rev",
                    "fineqpy_rev",
                    "wsekpwr_rev",
                    "weasoff_rev",
                    "wexashr_rev")


reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence")


model_1 <- lmer(eqwrkbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
              data = data, REML = FALSE,weights=pspwght, 
              control = lmerControl(optimizer = "optimx", 
                                    optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqpolbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(eqmgmbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(eqpaybg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Equality_family life",
                                            "Equality_politics",
                                            "Equality_work",
                                            "Equality_economy"),
          show.ci=F,show.se=T,collapse.se=T,file="h2a_1.doc")


model_1 <- lmer(eqparep_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqparlv_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(freinsw_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(fineqpy_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Policy_parliament",
                                            "Policy_parental leave",
                                            "Policy_insult",
                                            "Policy_payments"),
          show.ci=F,show.se=T,collapse.se=T,file="h2a_2.doc")


model_1 <- lmer(wsekpwr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(weasoff_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(wexashr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Stereotype_control",
                                            "Stereotype_offend",
                                            "Stereotype_exaggerate"),
          show.ci=F,show.se=T,collapse.se=T,file="h2a_3.doc")

##H2b

dependent_vars <- c("eqwrkbg",
                    "eqpolbg",
                    "eqmgmbg",
                    "eqpaybg",
                    "eqparep_rev",
                    "eqparlv_rev",
                    "freinsw_rev",
                    "fineqpy_rev",
                    "wsekpwr_rev",
                    "weasoff_rev",
                    "wexashr_rev")


reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence","R's gender X I's gender")


model_1 <- lmer(eqwrkbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqpolbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(eqmgmbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(eqpaybg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Equality_family life",
                                            "Equality_politics",
                                            "Equality_work",
                                            "Equality_economy"),
          show.ci=F,show.se=T,collapse.se=T,file="h2b_1.doc")


model_1 <- lmer(eqparep_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqparlv_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(freinsw_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(fineqpy_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Policy_parliament",
                                            "Policy_parental leave",
                                            "Policy_insult",
                                            "Policy_payments"),
          show.ci=F,show.se=T,collapse.se=T,file="h2b_2.doc")


model_1 <- lmer(wsekpwr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(weasoff_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(wexashr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender*gender_int+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Stereotype_control",
                                            "Stereotype_offend",
                                            "Stereotype_exaggerate"),
          show.ci=F,show.se=T,collapse.se=T,file="h2b_3.doc")


###H3

dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")

reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence","GII","I's gender X GII")

model_1<- glmer(unfair_med~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gii_sc + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_2<- glmer(unfair_work~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gii_sc + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_3<- glmer(unfair_police~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gii_sc+ (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_4<- glmer(unfair~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*gii_sc + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))


tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Medical treatment",
                                            "Work",
                                            "Police encounter",
                                            "Any domain"),
          collapse.ci=T,file="h3_1.doc")
###cont

reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence","GII","I's gender X GII")


model_1 <- lmer(eqwrkbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqpolbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(eqmgmbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(eqpaybg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Equality_family life",
                                            "Equality_politics",
                                            "Equality_work",
                                            "Equality_economy"),
          show.ci=F,show.se=T,collapse.se=T,file="h3_2.doc")


model_1 <- lmer(eqparep_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqparlv_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+
                  intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(freinsw_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(fineqpy_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Policy_parliament",
                                            "Policy_parental leave",
                                            "Policy_insult",
                                            "Policy_payments"),
          show.ci=F,show.se=T,collapse.se=T,file="h3_3.doc")


model_1 <- lmer(wsekpwr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(weasoff_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(wexashr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*gii_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Stereotype_control",
                                            "Stereotype_offend",
                                            "Stereotype_exaggerate"),
          show.ci=F,show.se=T,collapse.se=T,file="h3_4.doc")

###seats


dependent_vars <- c("unfair_med","unfair_work","unfair_police","unfair")

reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence","Parliamentary seats","I's gender X Seats")

model_1<- glmer(unfair_med~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*seats_sc + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_2<- glmer(unfair_work~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*seats_sc + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_3<- glmer(unfair_police~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*seats_sc+ (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

model_4<- glmer(unfair~ gender + age_sc + edulvlb_cont_sc + domicil + marital +
                  child + work + hincfel + relig_reverse_sc + lrscale_sc + 
                  gender_int + intage_sc +
                  vdcond + prewhp_new + gender_int*seats_sc + (1|intnum_new) + (1|psu) + (1|cntry),
                data = data, family = binomial, weights = pspwght,
                control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))


tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Medical treatment",
                                            "Work",
                                            "Police encounter",
                                            "Any domain"),
          collapse.ci=T,file="h3_b1.doc")


###cont

reg_label <- c("Intercept","Respondent's gender (ref.: male)","Respondent's age",
               "Education","Settlement size",
               "Married","Lives/lived with children","Ever done payed work","Income",
               "Religious attendance","Left-right position",
               "Interviewer's gender (ref.: male)",
               "Interviewer's age","Interview outside (ref.: at home)","Video intreview",
               "Spousal or partner presence","Parliamentary seats","I's gender X Seats")


model_1 <- lmer(eqwrkbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqpolbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(eqmgmbg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(eqpaybg~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Equality_family life",
                                            "Equality_politics",
                                            "Equality_work",
                                            "Equality_economy"),
          show.ci=F,show.se=T,collapse.se=T,file="h3_b2.doc")


model_1 <- lmer(eqparep_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(eqparlv_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+
                  intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(freinsw_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_4 <- lmer(fineqpy_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Policy_parliament",
                                            "Policy_parental leave",
                                            "Policy_insult",
                                            "Policy_payments"),
          show.ci=F,show.se=T,collapse.se=T,file="h3_b3.doc")


model_1 <- lmer(wsekpwr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_2 <- lmer(weasoff_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))


model_3 <- lmer(wexashr_rev~ gender+age_sc+edulvlb_cont_sc+domicil+marital+
                  child+work+hincfel+relig_reverse_sc+lrscale_sc+gender_int+intage_sc+
                  vdcond+ prewhp_new+gender_int*seats_sc+(1|intnum_new) + (1|psu)+ (1|cntry),
                data = data, REML = FALSE,weights=pspwght, 
                control = lmerControl(optimizer = "optimx", 
                                      optCtrl = list(method = 'nlminb')))

tab_model(model_1,model_2,model_3,model_4,
          pred.labels=reg_label,dv.labels=c("Stereotype_control",
                                            "Stereotype_offend",
                                            "Stereotype_exaggerate"),
          show.ci=F,show.se=T,collapse.se=T,file="h3_b4.doc")
